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Abstract 

Systems of linear equations with Toeplitz coefficient matrices arise in many important appli- 
cations. The classical Levinson algorithm computes solutions of Toeplitz systems with only 
0(n 2 ) arithmetic operations, as compared to 0(n 3 ) operations that are needed for solving 
general linear systems. However, the Levinson algorithm in its original form requires that all 
leading principal submatrices are nonsingular. In this paper, an extension of the Levinson 
algorithm to general Toeplitz systems is presented. The algorithm uses look-ahead to sk' ■> 
over exactly singular, as well as ill-conditioned leading submatrices, and, at the same time, 
it still fully exploits the Toeplitz structure. In our derivation of this algorithm, we make 
use of the intimate connection of Toeplitz matrices with formally biorthogonal polynomi- 
als. In particular, the occurrence of singular or ill-conditioned submatrices corresponds to 
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Space Research Association. 

f The research of this author was supported in part by Army contract number DAAL-03-90-G-0105 and 
in part by Cooperative Agreement NCC 2-387 between the National Aeronautics and Space Administration 
and the Universities Space Research Association. 
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a breakdown or near-breakdown in the standard recurrence relations for biorthogonal poly- 
nomials. We present new general recurrence relations that connect successive pairs in any 
given subsequence of all existing formally biorthogonal polynomials. These recurrences then 
immediately lead to the proposed look-ahead Levinson algorithm for solving Toeplitz sys- 
tems. Implementation details for this algorithm and operations counts are given. Numerical 
experiments for Toeplitz systems with ill-conditioned submatrices are reported. 


1 Introduction 


Matrices whose entries are equal along each diagonal are called Toeplitz matrices. In partic- 
ular, a general square Toeplitz matrix of order n + 1 is of the form 


T = It- 1 - 


to t-\ t - 2 * * * t _ n 

t\ to ' • « 

t 2 * • * * ' • : 



t-\ 

tl to _ 


( 1 . 1 ) 


where the entries t x are real or complex numbers. In this paper, we are concerned with the 
solution of systems of linear equations 


T x — b 

A n^n — u n 


( 1 . 2 ) 


with Toeplitz coefficient matrices (1.1). The task of solving Toeplitz systems (1.2) arises in 
many important applications, such as time-series analysis [13, 29], linear prediction [20, 36], 
spectral estimation [22, 24], system identification [26, 27, 30], Pade approximation [6, 18], 
and statistics [19]. 

There are classical fast algorithms for solving (1.2) that exploit the Toeplitz structure and 
require only 0(n 2 ) operations, as compared to 0(n 3 ) operations for general linear systems. 
These algorithms implicitly compute either an inverse triangular factorization of T n of the 
type 

VjT n U n = D n , (1.3) 

or a triangular factorization of T n of the form 

T n = Vj D n U n , (1.4) 

see, e.g., [12]. Here, in (1.3) and (1.4), t/„ and V n are unit upper triangular matrices, and D n 
is a diagonal matrix. The class of fast Toeplitz solvers based on (1.3) includes the Levinson 
algorithm [35] and its variants [13, 43, 48, 49]. Toeplitz solvers based on (1.4) are intimately 
connected with the classical work of Schur [39, 31], and they are called Schur-type methods. 
Algorithms in this second class were proposed by Bareiss [2], Rissanen [38], and others; we 
refer the reader to [31] and the references given there. 
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Note that, for a nonsingular matrix T n , triangular decompositions of the type (1.3) 
and (1.4) exist if, and only if, all leading principal submatrices of T n are nonsingular. Indeed, 
all classical fast Toeplitz solvers require that T n is strongly regular , i.e., the submatrices T m , 
m = 0, 1, . . . ,n — 1, are all nonsingular. In some, but by far not all applications that lead to 
Toeplitz systems the coefficient matrices T„ are Hermitian positive definite and thus guar- 
anteed to be strongly regular. However, Hermitian indefinite and non-Hermitian Toeplitz 
matrices are in general not strongly regular, and it cannot be excluded that singular or 
ill-conditioned submatrices occur. We remark that Hermitian indefinite Toeplitz systems 
arise, for instance, in spectral estimation [37], when inverse iteration is used to compute 
eigenvalues of Hermitian positive definite Toeplitz matrices, see [10, 22, 24]. 

It is well known that the Levinson algorithm and the Schur-type Toeplitz solvers can be 
extended to handle exactly singular leading principal submatrices, and numerous algorithms 
were proposed [11, 17, 20, 23, 40, 47], These algorithms are again based on triangular 
factorizations of the type (1.3) or (1.4), where now D n is a block diagonal matrix. More 
precisely, blocks of size h k > 1 in D n just correspond to h k - 1 consecutive singular leading 
submatrices of T n . 

In finite-precision arithmetic, it is not enough to skip only over exactly singular subma- 
trices, and a numerically robust Toeplitz solver also must be able to handle nonsingular, 
yet ill-conditioned leading principal submatrices. The literature on Toeplitz algorithms with 
this property is rather scarce. Sweet [41, 42] showed that, in principle, pivoting can be 
incorporated into the Bareiss algorithm, which allows to treat singular and nearly singular 
submatrices. However, there are some unresolved difficulties with this algorithm, such as the 
necessity for an a-priori choice of parameters, and the fact that a large pivot does not nec- 
essarily guarantee well-conditioned submatrices. Recently, Chan and Hansen [8] proposed a 
look-ahead modification of the Levinson algorithm for general Toeplitz systems. If a singu- 
lar or a nonsingular ill-conditioned submatrix occurs, then the algorithm looks ahead to the 
next well-conditioned leading submatrix, and instead of a standard Levinson step, a block 
step is performed. However, this look-ahead algorithm is not entirely satisfactory. The look- 
ahead strategy used in [8] requires condition number estimates for all leading submatrices, 
and this generates overhead of the order 0(n 2 ), even if it turns out that no block steps are 
necessary. Moreover, as we will demonstrate with an example in Section 7.2 below, there is 
a potential source for a breakdown of the algorithm if two or more consecutive block steps 
are performed. 

In this paper, we propose a look-ahead Levinson algorithm for general Toeplitz systems 
that is different from the one in [8]. In our derivation of this algorithm, we make use of the in- 
timate connection of Toeplitz matrices with formally biorthogonal polynomials (FBOPs). In 
particular, the occurrence of singular or ill-conditioned submatrices corresponds to a break- 
down or near-breakdown in the standard recurrence relations for biorthogonal polynomials. 
First, we derive new general recurrence relations that connect successive pairs in any given 
subsequence of all existing FBOPs. These recurrences then immediately lead to the proposed 
look-ahead Levinson algorithm for solving Toeplitz systems. 

The remainder of this paper is organized as follows. In Section 2, we introduce some 
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notation, and we give a formal definition of FBOPs associated with general bilinear forms. 
We then turn to bilinear forms induced by Toeplitz matrices, and in Section 3, we collect 
some basic properties of the corresponding FBOPs. In Section 4, we derive general recur- 
rence relations for FBOPs. In Section 5, we propose a look-ahead procedure for constructing 
FBOPs, and we describe some properties of this algorithm. In Section 6, we present our 
look-ahead Levinson algorithm for solving general Toeplitz systems. We give implementation 
details and operation counts, and we discuss the look-ahead strategy. In Section 7, we con- 
sider the look-ahead Levinson algorithm for the special case of Hermitian Toeplitz systems. 
Also, we show that the procedure proposed by Chan and Hansen has potential breakdowns. 
In Section 8, we report results of numerical experiments with Toeplitz matrices that have 
various kinds of ill-conditioned submatrices. Finally, in Section 9, we make some concluding 
remarks. 


2 Preliminaries 

In this section, we introduce some notation, and we give a formal definition of FBOPs 
associated with general bilinear forms. 


2.1 Notation 


Throughout the paper, all vectors and matrices are allowed to have real or complex entries. 
As usual, M T := [ ] , M := and M H := M denote the transpose, complex 

conjugate, and conjugate transpose, respectively, of a matrix M = [ ] . The vector norm 

||x|| := i/S is the Euclidean norm, and ||M|| := max||j:|| = i ||Mx|| is the corresponding 
matrix norm. For square matrices M £ C hxh , we use the following condition number: 


ri/|M|, ith = h 

l||M|H|M- 1 ||, if h > 1. 


( 2 . 1 ) 


Whenever we call a square matrix ill-conditioned , it is with respect to the condition num- 
b< (2.1). We denote by Ik £ the k x k identity matrix, by 


J k := 


0 

.1 0 


0 n 

1 o 

... o 


£ R 


kx k 


t .iekxk antidiagonal identity matrix, by Ojtxj the k x j zero matrix, and by 0* € R* the 
zero vector of length k. We will drop subscripts and simply write /, J, or 0 if the actual 
dimensions are apparent from the context. 

The set of all complex polynomials of degree at most n is denoted by 


V n := {</?(A) = <To + 0 iA + '" + <r n A n | 00,01, . . . ,u n £ C}, 
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and V is the set of all complex polynomials. We denote by d(p) the exact degree of p £ V, 
i.e., d{p) is the smallest integer n > 0 such that p £ V n . A polynomial p £ V n is called 
monic if it is of exact degree n with leading coefficient 1. For each p £ V , we define its 
reverse p by 

p(X) = A a(v) p(l/X). (2-2) 

Note that p is a polynomial of degree at most d(p). 

Capital Greek letters are always used to denote row vectors of polynomials in V, e.g., 

$ = [^0 H>\ ••• ] ■ ( 2 - 3 ) 

A vector of the type (2.3) is called a block of polynomials. The reverse $ of a block (2.3) is 
defined by 

$(A) = A”*4>(l/A), where n t := .max %,•). (2.4) 

Note that the entries of $ are again polynomials, and 

$ = [ A (n *~ 9(v,o)) ^o A (n * -9(vi)) v3i ••• A (n * -9(v ’> ) )£ 7 - ] . (2.5) 

Furthermore, if j = 0 in (2.3), then (2.4) reduces to the usual reverse (2.2) of a single 
polynomial. 

Finally, for each n, we denote by 

A„ = [ 1 A A 2 ••• A"] (2.6) 

the block whose entries are the monomials A', i = 0,1,..., n. We remark that, with (2.6) 
and from (2.4), any block (2.3) and its reverse can be represented in the form 

1> = A and $ = A n JU, (2.7) 

respectively. Here U £ is a matrix whose ith column just contains the coeffi- 

cients of the polynomial tpi, i = 0, 1, . . . , j. In particular, each p £ V n can be written in the 
form 

tp — A n ii for some u £ C n+1 . (2-8) 

2.2 Bilinear Forms and FBOPs 

A complex-valued functional 

(2.9) 

is called a bilinear form if 

{ipiCTiPi + C 2 P 2 ) = for V>> ¥>ii ^2 € V, 0\, <t 2 £ C, ^ 

+ T 2 tp 2 ,p) = + T 2 (tp 2 ,p) for all V> 1 , 02, V € V, T X , T 2 £ C. 

We stress that, in general, a bilinear form is not an inner product. Indeed, it is possible that 
a nonzero polynomials p has “norm” (p, p) = 0 or (p, p) < 0. 
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Nevertheless, it turns out to be useful to study polynomials that are orthogonal with 
respect to a given bilinear form (2.9). Next we give a definition of these formally biorthogonal 
polynomials (FBOPs). 

Definition 2.1 A monic polynomial p n € V n is called a right FBOP (with respect to the 
bilinear form (2.9)) of degree n if 

(if,<p n ) = 0 for all ifeVn-i- (2-11) 

A monic polynomial rp n € V n is called a left FBOP (with respect to the bilinear form (2.9)) 
of degree n if 

(rf n ,<p)=0 for all ip e V n - a . (2-12) 

A right or left FBOP p> n or ip n is said to be regular if it is uniquely determined by (2.11) 
or (2.12), respectively. 

REMARK. In general, regular FBOPs need not exist for every degree n; for instance, see 
Lemma 3.1 below. 

REMARK. Biorthogonal polynomials have been studied in various settings; we refer the 
reader to [3, 6, 32, 33, 46] and the papers cited therein. The notation “formally biorthogonal 
polynomials” goes back at least to van Rossum [46]. However, almost all the literature is 
concerned with cases where regular FBOPs of every degree n are guaranteed to exist or are 
assumed to exist. 

In the sequel, it will be convenient to use the following extension of the bilinear form (2.9) 
to blocks of polynomials of the type (2.3). More precisely, for any blocks 

$ = [<Po <pi ••• <Pj ] and = [if> 0 t/»i ••• i>k ) , 


we define 


<*,*> 


' (V>o, V>o) 
-(^fc.9 o) 


(V’o, <Pj) 


£ £(*+l)x(j+l) 


(V>Jfe, Vi) -1 


3 FBOPs Associated with Toeplitz Matrices 

Let {t,}^_ 00 be a given biinfinite sequence 1 of real or complex numbers, and let 

T n ~ [ U-j ], l j=o,...,n > n = 0, 1, < • • , 

be the associated family of Toeplitz matrices (1.1). The sequence {T n }£L 0 induces a bilinear 
form (xf,ip) on V x V as follows. Using the representation (2.8), for any two polynomials 

ip = A n u, if = A n v, u, v € C n+1 , 

Bn the case that only a finite sequence . . . , t n is given, we can always extend it to a biinfinite 

one, by simply setting ti := <_< := 0 for all t > n. 
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of degree at most n, n = 0, 1, . . . , we set 

<0,¥>) :=v T T n u . (3.1) 

From now on, we assume that {*,*) is the bilinear form defined by (3.1). Furthermore, the 
term FBOP always refers to formally biorthogonal polynomials with respect to this particular 
bilinear form. 

Next, we list some properties of {•,•). For the monomials v?(A) = X : and V’(A) = A’, we 
obtain 

(AS A = *, j = 0,1,..., (3.2) 

and hence the elements of {ti}°l_ oc are just the moments associated with (•,•). From (3.2) 
and the bilinearity relations (2.10), it readily follows that 

(Ai/>, Xp) = (rp,<p) for all (3.3) 

Furthermore, for all <p £ V and all j, k — 0, 1, ... , with d(p) + k — j > 0, we have 

(X k p,X j ) = {X dM+k ~ j ,p), ( 3 - 4 ) 

(AS A k (p) = {^\ d ^ +k - 3 ). (3.5) 

For example, to verify (3.4), we set n := d(p), and we represent p, p in the form 

y?(A) = X^o-.A’, £(A) = ^o-.A" - ', with a 0 ,<7j, . . . ,cr n € C. 

«=o »=o 

Then, with (2.10) and (3.2), we obtain 

(X k p, A J ) = £>*(*+„- i)-j = £>*(»+*->)-< = (A n+ * - S </>)• 

t =0 1=0 


Similarly, one can show (3.5). 

In view of (2.11), (3.1), and (1.1), a polynomial 

^ n (A) = A" + £u in A‘ 

1=0 

is a right FBOP of degree n if, and only if, its coefficients tz m satisfy 


Wo n 


' t-n ‘ 

U\n 


t-(n-l) 



• 

- Un—l,n - 


. t - 1 . 


(3.6) 


Analogously, by (2.12), (3.1), and (1.1), a polynomial 

V’n(A) = A" + £ VinX* 

i=0 
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is a left FBOP of degree n if, and only if, its coefficients u, n fulfill 

[ ^On ^ln • • • ^n— l,n ] T n — 1 = [ ^n— 1 • • • ^1 ] • (3-7) 

As an immediate consequence of (3.6) and (3.7), we have the following result. 

Lemma 3.1 The following conditions are equivalent: 

(i) A regular right FBOP <p n of degree n exists. 

(n) A regular left FBOP ip n of degree n exists. 

(in) The matrix T n _ i is nonsingular. 

4 Recurrence Relations for FBOPs 

In general, regular FBOPs <p n and 0 n need not exist for every n. We denote by 

{nj}^_ 0 , where 0 =: n 0 < nj < • • • < nj < ■ ■ ■ , (4.1) 

the sequence of all integers n for which regular FBOPs of degree n exist. Here, either J = oo 
or, if there are only finitely many regular FBOPs, J is an integer. We remark that, for n = 0, 
the conditions (3.6) and (3.7) are void, and thus y>o(A) = 1, V’o(A) = 1 are regular FBOPs 
of degree 0. Hence no = 0 is always included in (4.1). Note that, in view of Lemma 3.1, the 
sequence {n_, }j =1 just consists of all integers n > 1 for which T n _i is nonsingular. 

In this section, we derive recurrence relations for generating regular FBOPs corresponding 
to arbitrary subsequences of (4.1). 

4.1 The Classical Szego Recursions 

If nj = j and thus regular FBOPs of degree n exist for every n, then they can be generated 
by means of the celebrated Szego recursions [19, 3, 32]. 

Algorithm 4.1 (Szego recursions) 

0) Set <p 0 = V'o = 1, = (1,1). 

For n = 0, 1, . . . , do: 

1) Compute p n = (1,A<£>„), r n = (Ar/> n ,l). 


2) Set 




On = Pn/6n, <fn+l = A(^„ - V>nO„, 

fln = T n/&m V’n+l = A^n tfnPn- 

(4.2) 

3) Set 


<$n+l = Ml “ On^n)- 

(4.3) 
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We remark that the S n 's in (4.3) satisfy 

Sn = (V>n,¥>n). ( 4 - 4 ) 

The Szego recursions are no longer valid if regular FBOPs do not exist for every n. There 
are extensions of the relations (4.2) that connect consecutive pairs (f nj , t/’n, and 1 •> V’nj+i 
in the sequence (4.1) of all existing regular FBOPs. For example, such recurrences are given 
in [23] and, for the special case of Hermitian Toeplitz matrices {T n }£L 0 , in [11]. These 
recursions immediately lead to an extension of the Levinson algorithm that can skip over 
exactly singular leading submatrices. For the derivation of a more robust Levinson algorithm 
that can also skip over nonsingular, but ill-conditioned submatrices, we need more general 
recurrence relations that connect consecutive pairs in a subsequence of all existing regular 
FBOPs. In the next section, we present such general recursions for FBOPs. 


4.2 General Recursions for FBOPs 

Let {n it }£l 0 C {n_,}/ =0 be an arbitrary, but fixed, subsequence of (4.1). Here either K = oo 
or K is an integer. For simplicity, we set n k := n Jfc . Moreover, in view of (4.1), we can 
assume that, without loss of generality, n 0 = 0 is included in the subsequence. Therefore, 
we always have 

0 =: n 0 < ni < • • • < n* < • • • . (4.5) 

For all 0 < k < K, we set 

h k ■ ~ n,k + 1 n k and F k . | . . -> 

[ {n \ n k < n < n k+ i j, 

If /S' < oo, then we set nx+i •— oo and X# := {n | n > nx}- 

The goal then is to give recurrences for generating the regular FBOPs 

{^n*}fcLo and {VsjLo ( 4J ) 


if h k = 1, 
if h k > 1. 


(4.6) 


corresponding to the prescribed indices (4.5). To this end, we also construct additional monic 
polynomials 

(fin, V»» € Vn for all n € I k and all k. (4.8) 

The polynomials (4.8) are called inner polynomials. Note that the regular FBOPs (4.7) 
together with the inner polynomials (4.8) build two sequences of monic polynomials (y n } n= o 
and {rp n }%Lo that both span V. Of course, we still need to specify how to actually choose 
the inner polynomials. In order to obtain recurrence relations that involve as few as possible 
previous polynomials, it is crucial to construct the inner polynomials <^> n , n € X*,, as 
(juasi-FBOPs, in the sense that they satisfy the relaxed biorthogonality relations 

(V>,<^n)=0 for all tl> € V nk -\, 

(ip n ,<p) = 0 for all <f>eV nk -\. 


(4.9) 
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Next we introduce some further notation. For all 0 < k < K, we define blocks of 
polynomials 

:= [l^n* i Pn k + 1 i V , n*+ 1 *■’ V’n Ji+1 -l]- (4-10) 

If K < oo, then we also define infinite row vectors of polynomials 

$ (/0 :=[V?n*- Vn K +l •••], , I f(A) := [V>n K V'njf+l ••']• 

Moreover, for all 0 < k < K, we set 

A w :=[A n * A n * +1 ••• A n *+ 1-1 ] , (4.11) 

F (k) := (A (fc) ,<J> (fc) ), G w := (¥ k \A {k) ). (4.12) 

We remark that, with (4.11), the biorthogonality conditions (2.11) and (2.12) for the reg- 
ular FBOPs (4.7) and the relaxed biorthogonality relations (4.9) for quasi-FBOPs can be 
summarized as follows: the polynomials p n , tp n , n — 0, 1, . . ., are required to satisfy 

{A^ m \ <£>„) = 0 and (^ n ,A^ m ^) = 0 for all m = 0, 1, . . . , k(n) — 1, (4-13) 

where k(n ) is the integer such that n^( n ) < n < n*( n ) + i. 

Note that F W and defined in (4.12) are hk x hk matrices, with hk from (4.6). We 

will need the fact that these matrices are nonsingular. 

Lemma 4.2 If the inner polynomials (4.8) are constructed as gwasi-FBOPs, then the ma- 
trices F ( k ) and G ^ are nonsingular for all 0 < k < K . 

Proof. Suppose that F^ is singular. Then there exists a vector z £ C hk such that 

F^z = 0 and z 0, (4-14) 

and we set 

V - (4-15) 

Clearly, p is a monic polynomial of degree nk+ 1 , and since the polynomials in the block 
are linearly independent, we have <p ^ <p nk+1 - Using (4.13)-(4.15), we deduce that 

. , , , , . ,,, f 0, if 0 < m < k. 

= (AH, ^ ) + (AW, *»), = { F(i)z _ 0 _ , ( m - t 

Hence <p is a regular right FBOP of degree n^ + i, and this contradicts the uniqueness of 
regular FBOPs. 

Similarly, one shows that G W is nonsingular. 0 

For the formulation of our recurrence relations, we will also need the following quantities. 
If hk > 1, we define 

/*:=[/*,-. 0 k ,_ 1 ]F«‘ ) [°Y 1 ]. »:=[/*,- 1 0 k ,. 1 ](G«‘>) 7 '[°‘‘ 1 -']. (4.16) 
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and, if hk = 1, we set fk = 9k = 0- We remark that fk and gk just consist of the first hk — 1 
elements of the last columns of F and (G' ( * ) ) T , respectively. Thus we have 


p( k ) _ 


* fk 

* * 


and (G (fc) ) r = 



(4.17) 


where the elements * in the lower right corners in (4.17) are lxl. 

After these preparations, we can now state our recurrence relations for generating the 
regular FBOPs (4.7) corresponding to the prescribed indices (4.5), together with inner 
polynomials (4.8) satisfying the relaxed biorthogonality conditions (4.9). First, we set 
<£>_! = = A -1 and $d _1) = tyf -1 ) = 0. Then, for all —1 < k < K and n fc+1 < n + 1 < n* +2 , 

we set: 


<p n+1 = A <p n - ^ {k) a. 


* £ 
k *=n *+1 


if n + 1 = n k + 1 , 
if n + 1 € Jjt+i, 


(4.18) 


and 


where 


Otn 


= (GW)- t 



Pn — (1 ? 


= (F (k) r' 





Vn+l = - $“>/?„ - 


£ d" 1 *, 

t ,=n *+i 


if n + 1 = n fc+1 , 
if n + 1 £ Tfc+i , 


(4.19) 


(4.20) 


where 


fin = ( i ^)" 1 


■0*-i 
. T n 


^"n 


(AV’n,!), 


*n = (G^))^ 


’ 0 ‘ 
.9k. ‘ 


(4.21) 


Here, in (4.18) and (4.20), £ C are coefficients that can be chosen arbitrarily. 

Note that, by Lemma 4.2, the inverse matrices in (4.19) and (4.21) all exist. Moreover, 
we remark that the recurrence (4.20) can be equivalently formulated in terms of the reverse 
polynomials that appear in (4.18). The resulting relation is as follows: 




f if n + 1 — n.fe+ 1 , 

£ ( 1 - n) A n+ 1 ~*V> t , ifn + l€lfc+i. 


=Kk + l 


Of course, we still need to verify that the recursions (4.18)-(4.21) indeed generate the 
regular FBOPs (4.7). 


Theorem 4.3 Let {v?„}£L 0 and {V>n}£L 0 be the se 9 uences °f polynomials defined by the 
recurrence relations (4. 18)— (4.21 ) . Then, these polynomials satisfy the biorthogonality rela- 
tions (4.13). In particular, the polynomials }jt=o and {VOLo are the uniquely defined 
regular FBOPs corresponding to the prescribed indices (4.5). 
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Proof. We show (4.13) by induction on n. For n = 0, by (4.5), A:(0) = 0, and the condi- 
tions (4.13) are void. 

Now let n > 0, and assume that (4.13) holds for all polynomials ip 0 , ■ ■ ■ , '■Pn and 
ipo, il>i, . . . ,ip n - We need to show that <p n+i and V'n+i satisfy 


(A',9»+.)=0 

for all j = 0,1,. 

• • - 1, 

(4.22) 

«Vh,A’) =0 

for all m = 0, 1, . 

• • i ^fc+i 1 , 

(4.23) 


respectively. Here k is the integer defined by n k + 1 < n + 1 < n^+ 2 • To simplify notation, 
we set h' := /j*. — 1 and n' := n k +i — 1. Moreover, in the following, we always assume that 
j € {0,1,..., n'}. 

First, we consider (4.22). Writing in the form (2.5) and using (3.5), one readily 
verifies that 

= {¥ k \\ n ‘~ i ) T . (4.24) 

From (4.24) and (4.13), it follows that 

$ (fc *) = 0 for all j with h'<j<n'. (4.25) 

Using notations introduced in (2.6), (4.11), and (4.12), we can summarize (4.24) for the 
remaining indices 0 < j < h' as follows: 

(A^, $ w ) = (¥ k) , A W J) T = J(G (fc) ) r . (4.26) 


Next, we note that, in view of (3.3) and (4.13), we have 


(A J ,A^ n ) = 


( 1 , \ip n ) — p n , 

(A J_1 , <^n) = 0, 

(A^ -1 , V? n ) = 0, 
l (A J_1 ,^n), 


if j = 0, 

if 1 5; j < n ' and n + 1 6 l k +\ , 
if 1 < j < «jt and n + 1 = njt+i, 
if rik < j < n‘ and n + 1 = «jt+i • 


(4.27) 


Using the vector f k defined by (4.16) and (4.12), we can summarize the relations (4.27) for 
n k < j < n' and n + 1 = njt+i as follows: 


(A (i) ,A^> 



if n + 1 = n fc+1 . 


(4.28) 


With these preparations, it now readily follows that y n+ i satisfies (4.22). With (4.18) 
and (4.13), we obtain 


(A' , ,<r’n+i) = (A J , Ayj n ) - (\ 3 ,ty (k) )a n - ^ 


if n + 1 = n k+u 
if n + 1 G X*;+ 1 . 


(4.29) 
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By (4.27), (4.25), and (4.13), all terms on the right-hand side of (4.29) vanish and hence 
(\ : ,<p n +i) = 0 for all j in the range h' < j < n k , if n + 1 = n k+ 1 , and for all j in the range 
ti < j < n', if n+1 € J„+i. For; in the range 0 < j < h!, we use (4.29), (4.13), (4.27), (4.26), 
and the definition of a n in (4.19) to deduce that 


{Ah' , ^n+1 ) 


Pi 

LOv J 


- J{G (k) ) T Ct n = 0 . 


For the case that n + 1 = n* + i, it remains to prove (4.22) for all j in the range n k < j < n' . 
Here we use (4.29), (4.28), (4.25), (4.12), and the definition of fi n in (4.19) to verify that 


(A<‘>, V „ +1 ) = 


0 

fk J 


(A W ,^K = 


0 

L/fcJ 


F w u n = 0. 


This concludes the proof of (4.22). 

To show (4.23) we proceed similarly. First, using (3.4) and (4.13), one verifies the relations 


(* {k \A h>) = (FW) T J, 

($W, A J , ) = 0 for all j with h' < j < n' . 

With (3.3) and (4.13), we obtain 

' (Al/>n, 1) = T n , if J = 0, 

(A0 n , \ 3 ) = i 0, if 1 < j < n' and n + 1 G Ik+ 1 , 

.0, if 1 ^ j < n k an( i n + 1 = ftfc+i 


(4.30) 


(4.31) 


(AV«,A'*>) = [0 si], if n + 1 = n, +i . 


(4.32) 


Here is the vector defined by (4.16) and (4.12). From (4.20) and (4.13), we have 


«’„ + „V> = (A# n ,V) - /£($<*>, A>) - 


if n + 1 = n k+ 1 , 
if n + 1 € Tfc+i . 


(4.33) 


Finally, using (4.30)-(4.33) and the definitions of /?„ and u n in (4.21), one easily veri- 
fies (4.23). □ 

REMARK. Our derivation of general recurrence relations is different from the one used 
in [23, 11] to obtain special recursions for all existing FBOPs. The approach in [23, 11] is 
based on Iohvidov’s results [25] on the structure of exactly singular Toeplitz matrices, and 
it cannot be directly used to prove general recurrences of the type (4.18)-(4.21). 


5 An A^orithm for Constructing FBOPs 

In this section, we propose an algorithm for constructing regular FBOP based on the general 
recurrence relations derived in Section 4.2, and we describe some properties of this algorithm. 
We use the notation introduced in Section 4.2. 
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5.1 The Algorithm 

The algorithm generates two sequences and °f monic polynomials, where, 

for each n, c p n and xj) n are either regular FBOPs or inner quasi-FBOPs. As in (4.5), we 
use the indices n^, k = 0,1,..., to mark the regular FBOPs, and we always set n 0 = 0, 
ip 0 z= ip 0 = 1. The index n = 0, 1, . . . , is used as an iteration counter, where in the course of 
the nth iteration the algorithm generates the next pair of polynomials <f n+ \ and ifin+i- For 
each fixed n, we define l = l(n ) by 


n\ < n < n/+i. (5.1) 

Note that 1 = l(n) is just the number of the last pair of regular FBOPs </?„, and ip ni with 
degree < n. In addition to the blocks k = 0, 1, . . . , / — 1, defined in (4.10), we set 

$ (0 :=[¥>n, </>n,+i ••• ¥>»], ^ (,) :=[0 nj n,+l ••• i>n}- (5.2) 

We call the blocks and complete if n -f 1 = nj +1 ; in this case, the next polynomials 
<p n +i and ipn+i are constructed as regular FBOPs. If n + 1 < n/+i, then the blocks (5.2) are 
still incomplete ; in this case, the next polynomials <f n +i and V’n+i are constructed as inner 
quasi-FBOPs and added to tfdO and respectively. Finally, as in (4.11), (4.12), we set 

A (,) :=[A n ‘ A n,+1 ••• A n ] , F (,) :=( A (,) ,$ {,) ), G (l) := {¥ l) , A (,) ). (5.3) 

Using these notations, we can rewrite the recurrences (4.18)-(4.21) in the form of the fol- 
lowing algorithm. 


Algorithm 5.1 (Construction of FBOPs) 

0) Set y? 0 = 4>o = 1, $ (0) = * (0) = 1, F<°> = G<°> = (1,1), n 0 = 0, l = 0. 

For n = 0, 1, ... , do: 

1) Compute 

p n = (l,\<p n ) and T n = (AV>„,l). (5.4) 

2) Decide whether to construct tp n+ 1 and as regular FBOPs or as inner polynomials 
and go to 3) or 4), respectively. 

3) (Regular step) Set 


= (G*'>) 


-T 


0 1 

Pn 



Pn = ( F ( ' ) )- 1 




¥>n+l = A 

tfn+1 = Wn ~ * {,) Pn ~ * ( 'V, n • 


(5.5) 

(5.6) 


Set n ; + 1 = n + l, / = /+!, = 0, and go to 5). 
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4) (Inner step) Set 


* n = (G^)~ T 

0 n = (F ((-1) ) _1 


¥>»+i = A<^ n - ¥ l -Va n - £ 

«=nj 

V»n+1 = ~ £ C^V,- 


(5.7) 

(5.8) 


5) Set *W = [ *W ^ n+1 ] , *W = [ *W V>„+i ) , F(') = (AW, *W), GW = (*W, AW). 


Note that the description of Algorithm 5.1 is still incomplete, since no criteria for the decision 
in step 2) are given. We defer a discussion of this so-called look-ahead strategy to Section 6.3. 
Here, we only remark that, in view of Lemma 4.2, the polynomials v? n+ i and V’n+i can be 
constructed as regular FBOPs only if the following necessary condition holds true: 

F W and (jW are nonsingular if n + l=n; + i. (5-9) 

In particular, (5.9) guarantees that the inverse matrices in (5.5)-(5.8) all exist. 


5.2 Inverse Block Factorization of Toeplitz Matrices 

Next, we show that Algorithm 5.1 yields a factorization of T n of the form (1.3), where D n 
is now a block diagonal matrix. In the following, let n > 0 be arbitrary, but fixed, and let 
l = l(n) be the corresponding index defined by (5.1). 

Recall that, by Theorem 4.3, the polynomials 

Wi)U and {«"=„ (5.10) 

satisfy the conditions (4.13), which are equivalent to the following block biorthogonality 
relations: 

(vfd*^ <|>( m )) — 0 for all k ^ m, fc, m = 0, 1, . . . , /. (5.11) 

Moreover, we set 

/)(*>:=(¥<*>,$<*)), k = 0,1,...,/. (5.12) 

Here, in (5.11) and (5.12), tyW and are the blocks defined in (4.10), if k < /, and in (5.2), 
if k = l. In view of (2.8), the monic polynomials (5.10) can be represented in the form 

Vi = AjUj, where u 3 := [u 0 > uy ••• i,j 1 f € C J+1 , 

j . (5-13) 

xpj = AjVj, where Vj := [v 0 j t>i j ••• Uj-ip 1] € C J . 

Using (5.13) and the definition of (•, •) in (3.1), we can rewrite the relations (5.11) and (5.12) 
in the following matrix formulation: 


VjZMr, = D 


(5.14) 
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where 


U n :— [ Uij ]■ j =0 „ 


1 «oi * ' • u 0n 

0 1 : 

• • * ^7i— 1 ,n 

0 ••• 0 1 


and 


V n ], j_ 0 n 


1 voi 
0 1 

0 ••• 


«On 

^n— l,n 

0 1 


are unit upper triangular matrices, and 


O n :=diag(£><°\.D< 1 >,... > fl (l >) 


(5.15) 


(5.16) 


(5.17) 


is block diagonal. In other words, Algorithm 5.1 generates an inverse block factoriza- 
tion (5.14) of T n , where the columns of U n and V n are just the coefficients of the poly- 
nomials (5.10). Furthermore, the sizes of the blocks in (5.17) are given by 


D w e c hkXhk , where h k 


{ Tlk+l ~ ™k, 
n + l - ni, 


if 0 < k < l, 
if k = l. 


(5.18) 


In particular, h k = 1 for all k if only regular steps 3) are performed in Algorithm 5.1. In 
the sequel, the construction of a true block of size h k > 1 will be referred to as a look-ahead 
step, and h k will be called the length of the look-ahead step. 

Finally, we remark that, by (5.13), the blocks and have the following represen- 
tations: 


= A nfc+1 -i 

* {k) = A n * +1 -lV<*>, 


where U ^ := [u,y > 

where V <*> := [ v {j ]- =0 j=n*,...,n* +l -i > 


if 0 < k < l — 1, and 

$ (,) = An UV\ where U {1) := [u,j]. =0 n;>=B| 

♦WssAnVW, where V (,) := [v,j] i=0 >n;j=Tl| , 

if k = l. Note that, by (5.14), (5.17), and (5.20), we have 

D d) = (yd)fT n U ( °. 


(5.19) 


(5.20) 


(5.21) 


5.3 Further Properties 

In this subsection, we collect some further properties of Algorithm 5.1 that will be needed 
later on. 


Formally Biorthogonal Polynomials and a Levinson Algorithm 


17 


Using (5.20), (4.13), (5.3), and the definition of (■,•) in (3.1), one readily verifies that 

T.U » 1 = [°J$' ] , (V m ) T T n = [0 k , xn , <?<'> ] . (5.22) 

Next, we partition the matrices and V ^ in the form 

(7 (,) = [^!!] , U (,) = , where V (,) € C fc ‘ x \ (5.23) 

Note that and are unit upper triangular matrices, and in particular, they are 
nonsingular. With (5.21) and (5.23), it follows from (5.22) that 

F {,) = (V {l) )~ T D {l \ G (l) - D 0) {U (l) )~\ (5.24) 

By means of the relations (5.24), F^ and G^ can be obtained without evaluating the bilinear 
forms (A (i) , $ (,) ) and (^ (() ,A (,) ) in step 5) of Algorithm 5.1, provided that D (i) is available. 
It turns out that the elements of D W can easily be updated from step to step, see Lemma 5.2 
below. 

Recall that the recurrence coefficients £,- n ^ and in (5.7) and (5.8), respectively, are 
still arbitrary. From now on, we assume that they are chosen as 

£jn) = c (n) = 0 for all i,n. (5.25) 

Furthermore, it will be convenient to set 

4*> := (G«)- t f 0 *;-' ] . 4" := ( ] ■ (5.26) 

Note that, with (5.26), the coefficient vectors a n , /?„ in step 3) and step 4) of Algorithm 5.1 
are given by 

a n = p n o^\ Pn = T n 0i^ and a n = p n a\-i, Pn = T n Pll\, ( 5 - 27 ) 

respectively. 

By (5.12) and (5.2), the matrix is given by 

#« = [<*.*>,•)]«_ ^ 5 - 28 ) 

In the following lemma, we give formulas for updating the elements of (5.28) in the course 
of Algorithm 5.1. 

Lemma 5.2 Let <p n +i, V’n+i be the polynomials constructed in step n Algorithm 5.1. 
a) If <fn+ii V'n+i are constructed as regular FBOPs, then 


(ifrn+l , <p n +l) = (V’n, <fn) ~ Pn [ T„ 0 gf ] JV^CX^ - [r n 0 gj ]U (l) p n . 


(5.29) 
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b) If (p n +i , V’n+i are constructed as inner polynomials, then 

(V’m+l i V’n+i ) = 


and 


{^n+l i ) — 


^0,nj Pn ? 

m + 1 = ni, 


(5.30) 

(0m j ^n) VQ rn+lPn j 

if ni < m + 1 < 

n + 1, 


z/ m -f 1 = n/ ? 
i/ n/ < m + 1 < 

n + 1 . 

(5.31) 



Proof. First we show part a). Here the polynomials tf n +i and <Pn+\ are given by (5.6) 
and (5.5), respectively. Using (5.6), the biorthogonality conditions (4.13), (5.5), (3.3), and 
the first formula for a„ in (5.27), one readily verifies that 

(0n+l,9n+l) = (Wn ~ $ (<) /?» ~ V?„+l ) 

= = (AV>„,Av? n - ^ (,) Q !„ - $ (, Vn) 

= (V’n, V’n) - P„(AV>„, - (Xlfn, (5.32) 


In view of (4.31) and (4.32), we have 

(A0 n , A n ) = [ r n 0 gf). (5.33) 

Note that, by (5.20) and (2.7), = A n U^ and = A n JV^ l \ Together with (5.33), it 

follows that 


(AV>n,^ {0 ) = [r„ 0 gj ] JV^ l \ (A^ n , $<'>) = [t„ 0 gJ)U^. (5.34) 

By inserting (5.34) into (5.32), we obtain the relation (5.29). 

Now we turn to part b). Since the proofs of (5.30) and (5.31) are complete analogous, 
we will only show (5.30). We assume that y? n +i is constructed as an inner polynomial; note 
that ifn+i is given by (5.7), where, by (5.25), £,- n ^ = 0 for all i. Let n;<m+l<n + l. 
From (5.6), if m + 1 = ni, and (5.8) (with (j m ^ = 0 for all i), if m + 1 > n/, it follows that 
V> TO +i can be written in the form 

Vwi(A) = Vwi(0) + AV>m(A) + Axm(A) for some Xm € V ni - 2 - (5.35) 

Moreover, in view of the representation (5.13) of t/’m+i, we have 

VWi(O) = u 0 , m +i- (5.36) 

Using (5.7), (4.13), (5.35), (3.3), (5.36), and the formula of p n in (5.4), one easily verifies 
that 

(tf m+1 ,y n+ i) = {^m+i,A^ n - ^ (/_1) a n ) = (Vwi, A <p n ) 

= ^ Pn) H" (At/j m , \<p n ) + (A Xmi A (fi n ) 

= v 0,m+l^n (VVnj^n)* 


(5.37) 
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If m + 1 > n;, then (5.37) is just the desired relation (5.30). If m -f 1 = n/, then, by (4.13) 
and since n > n; > m, we have (xf mi^n) = 0, and thus (5.37) reduces to (5.30). 0 

Typically, Algorithm 5.1 will perform mostly regular steps 3), and then the formulas (5.5) 
and (5.6) simplify somewhat. Indeed, assume that Algorithm 5.1 performs regular steps in 
two consecutive iterations, i.e., 

n + 1 = ni + 1 = Hi + 1 and hi = n/+i - nj = 1. (5.38) 

Using (5.3), (5.2), (4.13), and (5.12), one easily shows that, in this case, 

F (‘) = G U ) = ^ ipn ' ) = D (‘)_ (5.39) 


With (5.39), (4.16), (5.5), and (5.6), we conclude that 


fl — — 0? gn — k'n — 0, Q n 


<0n, Vn) 


, Pn = 


(0n,<Pn)' 


(5.40) 


In particular, in recurrences for ip n +i and 0 n +i in (5-5) and (5.6) reduce to 

(fn+i = A <p n - 4> n oi n and V’n+i = A^ n - <p n fin, ( 5 -41) 

respectively. Furthermore, with (5.40), (5.26), (5.39), and since JV^ = Jv n = [1 ••• ] , 

it follows that the update formula (5.29) reduces to 


(V>n+l,V?n+l) = (0n, 



Pn^n \ 
(V’n.Vn) 2 / ' 


(5.42) 


Note that, in view of (4.4), the recursions (5.41) and (5.42) are identical to the update 
formulas (4.2) and (4.3). In other words, the nth step of Algorithm 5.1 just reduces to the 
nth step of the Szego Algorithm 4.1, if (5.38) is satisfied. In particular, for the case that 
n*+i = n k + 1 for all k = 0, 1, . . we have the following result. 


Corollary 5.3 If all polynomials are constructed as regular FBOPs, then Algorithm 5.1 
reduces to the Szego Algorithm 4.1. 

6 A Look-Ahead Toeplitz Systems Solver 

In this section, we present our look-ahead Levinson algorithm for solving general Toeplitz 
systems. We give implementation details and operation counts, and we describe the look- 
ahead strategy. 
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6.1 The Algorithm 

By means of (5.13), Algorithm 5.1 can be rewritten in terms of the coefficients vectors u n 
and v„ of the polynomials <p n and t/>„, n = 0, 1, . . . . Recall from Section 5.2 that these vectors 
are just the columns of the triangular matrices U n and V n in the factorization (5.14)-(5.17) 
of T n . We call u n and v n regular vectors if they are the coefficient vectors of regular FBOPs 
if n and V’m and we refer to u n and v n as inner vectors if they correspond to a pair of inner 
polynomials. 

In the following, we denote by s n , r n G C n+1 the vectors defined by the partitioning 

r.+i = f‘° *”1 (6.i) 

L^„ i n J 

of T n +i. In view of (5.13), we have 

A v?n = A n +i [ ° and A ip n = A n+1 ° . (6.2) 

u n J Vn 

From (6.1), (6.2), and the definition of {•,•) in (3.1), it follows that 

(1 , A ip n ) = s^Un and {AV>„, 1) = vjr n . (6.3) 

Also, note that (1,1) = to- Finally, using (5.13), (6.3), (5.19), (5.20), (2.7), (5.25), and the 
vectors a^\ defined in (5.26), we can rewrite Algorithm 5.1 as follows. 


Algorithm 6.1 (Inverse block factorization of general T n ) 

0) Set Uq — Vq — 1, l/(°) = VW = 1, F (0) = G(°) = t 0 , n 0 = 0, l = 0. 
For n = 0, 1, . . . , do : 

1) Compute 


Pn =slu n , T n = vjr n . 


2) Decide whether to construct u n +i and v n +i as regular vectors or as inner vectors 
and go to 3) or 4), respectively. 

3) (Regular step) Compute a} 1 *, p n , fi\ X \ v n by solving 

(GW = [°] , F'V =[",]. F'V=[;], (G<'>f„ n = [°] , (6.5) 


respectively , and set 


UU+1 ~ UnJ Pn 


JV® 

0 

JUW 


a) - 


o i \jw'n a(i) rv^n 
Vn+1 ~ Uni Tn l o J^‘ [ o J Un ' 

Set n«+i = n + 1, / = / + 1, U® = W = /?(') = G<0 = 0, and go to 5). 
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4) (Inner step) Set 



■ 0 ' 


rji/O-Di 

(1) __ 

■ 0 ■ 

^n+1 — 

.Un. 

Pn 

0 

&1- 1 ? u n*fl 

V n . 


T n 


JUV-'U 

0 


& 


(1) 

/ — 1 • 


5) Set 



^n+1 ? 



^n+1 i 


and update F^ l \ 


Note that Algorithm 6.1 only computes the inverse block factorization (5.14) of T n . Next, 
we discuss how to obtain solutions for nested Toeplitz systems 


T n x n — b n , ti 0,1,..* . 

Here the right-hand sides b n £ C n+1 are assumed to be nested, i.e., 


(6.7) 


b n +l — 


bn 

L ®n+l J 


, with <y n+ i € C, for all n. 


Recall that T n is guaranteed to be nonsingular for n = — 1, k = 1,2,...,/, and we only 

update the solution x n of (6.7) for these values of n. To this end, we simply need to insert 
the following procedure at the beginning of each regular step 3) in Algorithm 6.1: 


Set n' = ni — 1, and partition U^\ b n , and T n as follows: 






K = 


bn' 


T = 

) J n — 


T n . S 
L R T n -n'\ ’ 


( 6 . 8 ) 


where and a just contain the last n — n' rows of and b n , respectively. 

Compute y by solving 

F®y = cr - Rx n >, (6-9) 


and set 


x n = 


X n t 

lo n _„,J 


+ u^y. 


( 6 . 10 ) 


We need to show that x n given by (6.10) and (6.9) is the solution of (6.7). Indeed, by means 
of (6.8)-(6.10) and the first relation in (5.22), it follows that 


T n x n = 


Tn' 

R 


Xn' + T n U (l) y = 


bn' 

Rx n ' 


+ 


f(‘) 


(F {,) )- l (cr-Rx n ,) = bn 
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6.2 Implementation Details and Operation Counts 

We now discuss some implementation details for Algorithm 6.1, and we present operation 
counts. 

Except for contrived examples, Toeplitz matrices that arise in practical applications have 
at most a small number of consecutive ill-conditioned leading principal submatrices. Conse- 
quently, Algorithm 6.1 mostly performs regular steps. Typically, only a few true look-ahead 
steps occurs, and their length hi is usually small, mostly hi = 2. This justifies the follow- 
ing convention that we will use for the operation count: a computation that requires only 
arithmetic operations of order O(hf) or less is considered negligible. 

We now consider steps l)-5) of Algorithm 6.1 in more detail. Step 1) involves the 
computation of two inner products of vectors of length n + 1. It turns out that these are the 
only two inner products that are required during the nth iteration. This is exactly the same 
as in the classical Levinson algorithm for strongly regular matrices. 

The look-ahead strategy for the decision in step 2) will be described in Section 6.3. As 
we will see there, it only involves negligible work. 

Next we turn to step 3). Note that (G^) T and F (,) are hi x hi matrices, and, if hi > 1, 
we use Gaussian elimination to solve the four linear systems in (6.5). Recall from (4.17) that 
fi and gi are given as part of the last columns of F^ and (G^) T . If hi = 1, then, by (5.40), 

= v n = 0, and the two updates in (6.6) require two SAXPYs 2 with vectors of length n + 1. 
If hi > 1, then we first compute the two vectors 

JV m ct\ l) and JU {l) p\'\ (6.11) 

which costs 2 hi SAXPYs. The two updates in (6.6) then require 2(h/-f 1) additional SAXPYs. 

Step 4) always requires two SAXPYs. This is obvious if /q_i = 1. If hi - 1 > 1, we use the 
vectors (6.11) (with l replaced by / — 1), which were already computed in the course of the 
last regular step. 

In step 5), we need to update the matrices F (/) and G (J) . To this end, we first update 
Z)0) using Lemma 5.2, and then we compute F^ and G^ by means of (5.24). Note that, 
by (5.23), the triangular matrices f/ (,) and V in (5.24) are given as part of U (l) and 
Consequently, step 5) only involves negligible work. 

Finally, we turn to the procedure for updating solutions of Toeplitz systems (6.7). To 
compute the right-hand side of (6.9), we need to generate Rx' n , which involves h { = n - ri 
inner products. The computation of the vector x n in (6.10) requires another hi SAXPYs. 
Note that x n is only updated once within each cycle of h t steps. Thus, in the average, the 
update procedure requires one inner product and one SAXPY per nth step. 

In Table 1, we summarize the operation counts for one step of Algorithm 6.1, and for the 
updating procedure for solutions of general Toeplitz systems (6.7). 


2 A SAXPY operation is z = x + ay, where x and y are vectors and a is a scalar, see, e.g., [16]. 
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regular step 

inner step 

update of x n 


with hi = 1 

with hi > 1 



inner products/step 

2 

2 

2 

i 

SAXPYs/step 

2 


2 

i 


Table 1: Operation counts for Algorithm 6.1 

6.3 The Look-Ahead Strategy 

The look-ahead strategy is crucial both for the accuracy and efficiency of Algorithm 6.1. 
Its main purpose is to skip over ill-conditioned leading principal submatrices in order to 
avoid breakdowns and numerical instabilities. However, as is obvious from the operation 
counts in Table 1, it is more expensive to perform a look-ahead step of length hi > 1 than 
hi classical Levinson steps. Therefore, for the sake of efficiency, it is desirable to perform 
look-ahead steps only when necessary. The look-ahead strategy is implemented through 
the criteria that are used in step 2) of Algorithm 6.1 to decide whether the next vectors are 
constructed as regular or inner ones. There are two quantities, denoted by k(T;) and r) n , that 
are monitored throughout the algorithm. Both of them are obtained from local information 
only. In particular, we do not need to estimate the condition number of the current leading 
principal submatrix T n . The decision about building regular or inner vectors is then based 
on a comparison of k(T^) and Tj n with two threshold parameters CDND (> 0) and GFACTOR 
(> 1), respectively. The algorithm dynamically determines COND and GFACTOR, and the only 
input that is required from the user is the number h m ax , which is the maximal length of a 
look-ahead step the algorithm is allowed to perform. Note that, in view of (5.18), we then 
have 

hk < h m&x for all k, 

and Algorithm 6.1 generates an inverse block factorization (5.14) with blocks of size < /i max . 
We remark that Algorithm 6.1 reduces to the classical Levinson algorithm if we set h max = 1. 

The initialization phase of the algorithm is as follows. As the first block, we set D = 
T m , where T m is the matrix with the smallest condition number 3 /c(T m ) among Th, h = 
0, 1, . . . , h mAX - 1. Then we build the next vector u m+ i and v m+ i as regular vectors and we 
set rii = m + 1. Furthermore, we initialize COND to be this smallest condition number /c(T m ), 
and we set GFACTOR to 1. 

Now we consider a general iteration step of Algorithm 6.1. As we mentioned before, for 
the sake of efficiency, the look-ahead Toeplitz solver should build as many regular vectors as 
possible. Therefore, in each iteration step, we first pretend that u n+ i and v n+i can actually 
be constructed as regular vectors. Recall from (5.9) that, for a regular step, it is necessary 
that F W and are nonsingular. To check this condition, we compute the matrix 

r (0 := 


3 Recall the definition of of k in (2.1). 
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where and V ^ are the triangular matrices given by (5.23). Note that, by (5.24), we 
have 

r(o = F (,) (i/ (0 ) -1 = (v {1) )~ t g {, \ 

and hence F M and G ^ are nonsingular if, and only if, is nonsingular. We now check 
whether 

K(r (,) ) < 2*C0ND. (6.12) 

If (6.12) is not satisfied, then we go to step 4) in Algorithm 6.1, and we build u n+1 and u n+1 
as inner vectors. To justify the criterion (6.12), recall that, in view of Lemma 3.1, T n is 
required to be nonsingular for a regular step. Actually, at the end of this section, we will 
point out that /c(r^) is closely related to n(T n ). 

If (6.12) is satisfied, then we compute the quantities a n , //„, /?„, and i/ n , using (6.5) and 
the first two relations in (5.27), and we set 

T) n : = max{||a n ||i,||^„||i,||/3 n ||i,||i/„||i}. 

Here ||x||i := |6| H f- |£fe| denotes the 1-norm of a vector x = [£i • • • (h ] T € C h . Then, 

we check whether 

r) n < 2*GFACT0R. (6.13) 

If the criterion (6.13) is not satisfied, we proceed with step 4) and construct u n+1 and t> n +i 

as inner vectors. The justification for (6.13) is as follows. Note that, by (5.14), 

T„ +1 = V-^D^U^, and T^ = U-^Dl^V-y (6.14) 

If one would compute the decompositions (6.14) of T n + 1 or T^ +x directly by means of Gaussian 
elimination, then pivoting would be used to ensure that the size of the off-diagonal elements 
of U~h and V~+ x is bounded by 1. Indeed, this is the key to numerical stability of Gaussian 
elimination. Recall that the elements of the strictly upper triangular parts of U~l x and 
are just the multipliers in Gaussian elimination. Now, for Toeplitz matrices pivoting would 
destroy the Toeplitz structure. Roughly speaking, the look-ahead Algorithm 6.1 performs 
a true look-ahead step of length hi > 1, whenever one would encounter a small pivot in 
Gaussian elimination. Consequently, the look-ahead strategy should also guarantee that 
off-diagonal elements of U~l x and are not too large. Since 

U~l x = I + C + C 2 + hC rn+1 , where C:=I-U n+ 1 , 

a large off-diagonal element of f/ n+ i usually leads to a large off-diagonal element of U ~+ A 
similar conclusion also holds for V ~+ x . Therefore, in each step of the Algorithm 6.1, we limit 
the growth in the newly generated off-diagonal elements of f / n + i and V n + X , which are just 
the components of u n+ i and v n+1 . This is the purpose of imposing check (6.13). 

If both criteria (6.12) and (6.13) are satisfied, then we proceed with step 3) in Algo- 
rithm 6.1 and construct u n+J and u„ + i as regular vectors. 

It cannot be excluded that the algorithm has reached the maximal block size /i max , but 
the two checks for building the next vectors as regular ones are still not satisfied. More 
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precisely, the algorithm has built a pair of blocks and of size h m ax , both starting 
with index n; and ending with index n; + /i max — 1, and the criteria (6.12) and (6.13) for 
constructing u n ,+^ miX and as regular vectors are not fulfilled. In this case, we adjust 

the values of the threshold parameters COND and GFACTOR in such a way that the criteria for 
building regular vectors are satisfied within the maximal look-ahead size h max . To this end, 
we first determine a step size h € {2, 3, . . . , h m&x } such that 

<c(if) + irS'W-. (s-is) 

is minimal. Here F { '\ j = denotes the j x j leading principal submatrix of 

r('). We then set COND to be the condition number of the corresponding T \ and the 
second parameter GFACTOR is set to the value of the corresponding r) nt +h-i- This choice of 
h guarantees that the vectors with index n/ + h can be constructed as regular vectors. The 
motivation for the choice (6.15) is as follows. From the above discussion, it is clear that the 
goal is to minimize simultaneously k(T^) and T] ni +h-i] this is exactly what (6.15) attempts 
to ensure. The weighting factor |r^| in (6.15) was chosen based on extensive numerical 
tests, and the choice (6.15) was found to work satisfactory in practice. With the described 
look-ahead strategy, the algorithm can be expected to produce accurate solutions of Toeplitz 
systems as long as the coefficient matrix has at most ft max - 1 consecutive ill-conditioned 
leading principal submatrices. 

We conclude this section with a discussion of the connection of the condition numbers of 
r6) and T n . We set n' := n; — 1, and we partition T n as follows: 


T = 

J. n — 



S' 

t. ‘ 


Here T n > is the last well-conditioned Toeplitz submatrix. With this notation, T (,) is, in fact, 
the Schur complement of T n < in T n , and we have the decomposition 


■ I 

O' 

T n , 

0 ' 

\I 

T„-'S1 

[RTJ 

/. 

. 0 


Lo 

/ . 


which implies that 


r<'> = [o 


/; 


-RT-, 





-T^Sl [O' 

I J [/. 


= [-RT~, 


-i 



(6.16) 


From (6.16), it follows that 





• Iir„||. 
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It can be verified that 


-RT-.' /hi < yrrpHJ^TPvfF, 

/ 


< v / > + (Il5||<c(r n .)/||r„. 


Therefore, setting 

e ■■= \A + (||B||K(7;.)/l|r„.||)Vi +(lisil*(r».)/lir„.||) 2 . 


we have 

l|r (0 || < p||T n ||. 


(6.17) 

Similarly, we can show 

|(r ( '>)-'| < e |r-'| 

. 

(6.18) 

To get bounds for the condition numbers, we distinguish two cases, 
is of size > 1. Then, by combining (6.17) and (6.18), we arrive at 

First, assume that P*' 


*(r.) < o J -c(r<"). 

With the same technique, we can show that 

ie(r<'>) < Q 2 K(T n ), 


and hence 

K(T n )/e 2 < k(T (,) ) < Q 2 K{T n ). (6.19) 

For the case that is a scalar, from (6.17) and (6.18), we obtain 

-^ n r<«(r (,) )<s||T- 1 ||. (6.20) 

Roughly speaking, the inequalities (6.19), respectively (6.20), state that, if T n > is well condi- 
tioned, then T n is well conditioned if, and only if, is well conditioned. 


7 The Special Case of Hermitian Toeplitz Matrices 

In this section, we consider the special case of Hermitian Toeplitz matrices. 

7.1 A Look-Ahead Levinson Algorithm 

Suppose that the elements of the biinfinite sequence {<,}“_«, satisfy f_, = t x for all i. Then 
the Toeplitz matrices (1.1) are all Hermitian, i.e., 

T n = T„, n = 0,1,.... (7.1) 
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Using (7.1) and the definition of (•,•), one easily verifies that regular FBOPs and inner 
quasi-FBOPs associated with the bilinear form {•,•) are connected by 

V’n = <?„- ( 7 - 2 ) 

Here ip denotes the polynomial whose coefficients are just the complex conjugates of the 
coefficients of <p. In other words, the coefficient vectors u n and v n of <p n and xp n satisfy 

v n = u n . (7.3) 

Consequently, for Hermitian Toeplitz matrices, Algorithm 6.1 simplifies, since we only need 
to update the vectors u n . Furthermore, note that, in (6.1), s n = r n, an d together with (7.3) 
and (6.4), we obtain 

r„ = v„r n = u„r n = {r„u n ) H = p n . (7.4) 

In view of (7.2), the matrices (4.12) are now connected by 

F (k) = (G {k) ) H . (7.5) 


Finally, using (7.3)-(7.5) and setting 


7 w „<>> = 


-1 


we obtain from Algorithm 6.1 the following look-ahead Levinson algorithm for Hermitian 
Toeplitz matrices. 


Algorithm 7.1 (Inverse block factorization of Hermitian T n ) 

0) Set uq = 1, UW = 1, F<°> = to, n 0 = 0, / = 0. 

For n = 0, 1, . . . , do: 

1) Compute r n = u„r n . 

2) Decide whether to construct u n+J as a regular vector or as an inner vector 
and go to 3) or 4), respectively. 

3) (Regular step) Compute 7 f 1 *, p n , by solving 





respectively , and set 



■ 0 ' 



fl) 


Un+l = 

,n„. 

- T n 

0 

h - 

. 0 J 


Pn- 


Set n ;+1 = n + 1, / = / + 1, U (,) = 0, and go to 5). 
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4) (Inner step) Set 


5) Set 


and update F“' 

7.2 A Counterexample 

Chan and Hansen proposed a look-ahead Levinson procedure that is different from our al- 
gorithm. Their method was first presented for the special case of real symmetric Toeplitz 
matrices [9], and then extended to general nonsymmetric Toeplitz systems in [8]. The deriva- 
tion of their algorithm is actually based on the assumption that only isolated look-ahead 
steps occur, which are preceded and followed by standard Levinson steps. However, in 
general it cannot be excluded that T n has two or more consecutive blocks of singular or 
ill-conditioned leading principal submatrices, which then requires two or more consecutive 
look-ahead steps. In both papers [9] and [8], this case is treated separately, and special 
recurrences are derived for handling two consecutive look-ahead steps, see [9, Theorem 3] 
and [8, Theorem2]. However, the proposed approach involves division by the first component 
Ci of a coefficient vector c, see [9, Equation (5.9)] and [8, Equation (50)]. Unfortunately, it 
is not guaranteed that Ci ^ 0, and thus division by 0 can occur. Indeed, we now present 
a real symmetric Toeplitz matrix for which Ci = 0. Consequently, the look-ahead Levinson 
algorithm proposed by Chan and Hansen can break down for general Toeplitz systems, as 
well as for the special case of Hermitian Toeplitz matrices. 

Consider the 7x7 Toeplitz matrix 

■0 1 0 1 1 oi- 
io 1 0 1 1 0 
0 10 10 11 

T 6 = 1 0 1 0 1 0 1 . (7.6) 

110 10 10 
0 110 10 1 
.10 110 10 . 

This matrix has exactly singular leading principal submatrices T 0 , T 2 , T 3 , and T 4 , while the 
remaining submatrices 7\, T 5 , and T 6 , are all nonsingular with condition numbers k(T\) = 1, 
/c(T s ) « 18.1, and k(T 6 ) « 7.2. In particular, T\ is optimally conditioned, and the algorithm 
in [9] starts with a look-ahead step of length 2, followed by another look-ahead step of 
length 4 or 5. For the update of quantities corresponding to the second look-ahead step, 
the algorithm [9] requires division by the first component c\ of a vector c that, using the 
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Figure 1: Condition number profile of a 64 x 64 Toeplitz matrix 
notations from [9], is given as follows: 

0 ' 

-1 
1 ' 

0 . 

Thus we have cj — 0, and the algorithm breaks down in the course of the second look-ahead 
step. 


( 2 ) _ 


Ta = 


r 

i 

0 

1 


J/2 = 


O' 

i 


Ra 


10 111 (2) r?T T 

[o i i op c= - rS ~ R < J y> = 


8 Numerical Experiments 


In this section, we report results of numerical experiments with the look-ahead Algorithms 6.1 
and 7.1, and the classical Levinson algorithm. All computations were carried out using 
MATLAB on a DEC 3100 workstation with machine precision of order 0(1O -16 ). For all 
the examples, we generated the right-hand side b n such that the vector of all l’s is the exact 
solution x exac t of T n x n = b n . In the following, we always list the relative error defined as 


relative error = 


1 1 -Z-compt 3- exact | 


(Inexact || 




where Xcompt is the computed solution. 

EXAMPLE 1. This test set consists of 100 nonsymmetric 64 x 64 matrices with at least 
one ill-conditioned leading principal submatrix. The off-diagonal entries t x of these matrices 
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Figure 2: Histogram of the relative errors for the classical Levinson algorithm, Example 1 


n + 1 

classical 

look-ahead 

overhead 

64 

200 

12096 

119400 




Table 2: Average number of multiplications for Examples 1 and 2 

were generated as random numbers in [—1,1], and to was then chosen so that at least on 
submatrix is ill-conditioned. A typical condition number profile of a matrix in this test set 
is shown in Figure 1. First, we use the classical Levinson algorithm, and in Figure 2, we 
show — in the form of a histogram — the relative errors for the 100 test matrices. We see that 
the relative errors are rather poor. In Figure 3, we plot the relative errors for the look-ahead 
Algorithm 6.1 (with h mtkX = 2). We see a substantial improvement of the relative errors. 

EXAMPLE 2. This test set consists of 100 nonsymmetric 200 x 200 matrices with at least one 
ill-conditioned leading principal submatrix. The matrices were generated as in Example 1. In 
Figure 4 and Figure 5, we plot the histograms of the relative errors for the classical Levinson 
algorithm and the look-ahead Algorithm 6.1 (with h mtLX = 2), respectively. In Table 2, we 
list, for both Example 1 and Example 2, the average number of multiplications required 
to solve one system in the test set by the classical Levinson algorithm and the look-ahead 
algorithm, and we state the corresponding overhead for the look-ahead algorithm. 

EXAMPLE 3. This test set consists of symmetric matrices given by 

n. 


t 0 = e, U = t-i = { l/2)‘, 


* = 1 , 2 ,..., 


(8.1) 
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Figure 3: Histogram of the relative errors for the look-ahead algorithm, Example 1 



Figure 4: Histogram of the relative errors for the classical Levinson algorithm, Example 2 





32 


Roland W. Freund and Hongyuan Zha 



Figure 5: Histogram of the relative errors for the look-ahead algorithm, Example 2 

These are special cases of a class of Toeplitz matrices called Kac-Murdock-Szego (KMS) 
matrices (with p = 1/2), see [28]. The eigenvalues of the KMS matrices (8.1) can be eas- 
ily computed [44]. It turns out that, for e = 0, each third principal submatrix, i.e., T 3m , 
m = 0,1,..., is exactly singular, while the remaining submatrices are well conditioned. 
Consequently, if t is set to a small number, then every third principal leading submatri- 
ces is ill-conditioned. The KMS matrices (8.1), with e = 10 -14 , were also used as test 
examples in [9], and we chose the same parameter t = 10 -14 . We ran the classical Levin- 
son algorithm and the look-ahead Levinson Algorithm 7.1 for KMS matrices T n of order 
n + 1 = 15, 30, 60, 120, 240, 480. The relative errors and the number of multiplications are 
listed in Table 3. We see that the look-ahead procedures yields solutions with nearly full 
accuracy. Considering that one third of the iteration steps are inner steps, the overhead of 
the look-ahead procedure is still reasonable. 

EXAMPLE 4. This test set consists of three small nonsymmetric matrices that have at least 
two consecutive ill-conditioned leading principal submatrices. To list the entires of 

a Toeplitz matrix (1.1), T n , of order n + 1, we use the following convention: 

t := to, t\ l— [f— 1 2 • • • t-n ] i • = [ ^1 ^2 ' ‘ ' In ] • 

EXAMPLE 4.1. This is a 5 x 5 matrix whose entries are listed in Table 4. As the condition 
number profile in Table 5 shows, this matrix has two consecutive ill-conditioned leading 
principal submatrices. 

EXAMPLE 4.2. This is a 6 x 6 matrix with three consecutive ill-conditioned leading principal 
submatrices. Its entries are listed in Table 6, and its condition number profile is shown in 
Table 7. 
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relative error 

multiplications 

n + 1 

classical 

look-ahead 

classical 

look-ahead 

overhead 

15 

0.0191 

1 .20“ 15 

630 

960 

52% 

30 

0.0184 

1.79 -15 

2610 

3870 

48% 

60 

0.0185 

1 .98“ 15 

10620 

15340 

46% 

120 

0.0186 

4.61 — 15 

42840 

62280 

45% 

240 

0.0186 

6.85' 15 

172080 

248333 

44% 

480 

0.0187 

3.69" 14 

689760 

995853 

44% 


Table 3: Test results for the KMS matrices, Example 3 


t 

tl 

t2 

-1.00000000000001 

1.27324683138786 

0.78539366864947 


-1.62115749363923 

3.41046741401696 


1.06413364195684 

-17.92422495778239 


1.21785304238395 

38.20692196916536 


Table 4: Elements of the matrix in Example 4.1 


m 

0 

1 

2 

3 

4 

K(T m ) 

1.00 

1 . 25e+15 

7.46e+14 

4.00e+l 

4.70e+2 


Table 5: Condition number profile of the matrix in Example 4.1 


t 

tl 

t2 

-0.99999999999998 

1.05288024249153 

0.94977563415339 


-1.10855680502906 

3.85673107101965 


1.16717755769466 

-13.61721591570147 


-2.22889818997626 

3.81850412563076 


4.51853189291597 

73.05176317918625 


Table 6: Elements of the matrix in Example 4.2 
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m 

0 

1 

2 

3 

4 

5 

<T m ) 

1.00 

7.70e+14 

6 . 33e+14 

3.91e+15 

4. 00e+l 

4 . 80e+2 


Table 7: Condition number profile of the matrix in Example 4.2 


t 

5 












tl 

-1 

6 

2 

5.697 

5.850 

3 

-5 

-2 

-7 

1 

10 

-15 

t2 

1 

-3 

12.755 

-19.656 

28.361 

-7 

-1 

2 

1 

-6 

1 

-0.5 


Table 8: Elements of the matrix in Example 4.3 

EXAMPLE 4.3. This example is taken from [41], and it was also used in [8]. It is a 13 x 13 
matrix with five consecutive ill-conditioned leading principal submatrices. Its entries are 
listed in Table 8, and the condition number profile is shown in Table 9. 

We ran the classical Levinson algorithm and the look-ahead Algorithm 6.1 for the three 
matrices in Example 4. The relative errors are presented in Table 10. 

EXAMPLE 5. Finally, we consider the symmetric 7x7 matrix (7.6) from Section 7.2. The 
classical Levinson algorithm breaks down for this matrix, while the look-ahead Algorithm 7.1 
generates the exact solution. Next, we perturb the matrix slightly by adding a Toeplitz 
matrix with random entries in [— 10 -14 , 10 -14 ] to (7.6). The condition number profile of 
the resulting matrix is shown in Table 11. Note that none of the submatrices are exactly 
singular. Nevertheless, the classical Levinson algorithm still breaks down. The look-ahead 
Algorithm 7.1 computes a solution with relative error 1.33e-14. 


9 Concluding Remarks 

We studied formally biorthogonal polynomials (FBOPs) for bilinear forms induced by general 
Toeplitz matrices. We presented new recurrence relations that connect successive pairs in 
any given subsequence of all existing FBOPs. Based on these recursions, we proposed a 


m 

0 

i 

2 

3 

4 

5 

K(T m ) 

2 . 00e-l 

i 

1.63 

9 .45e+5 

2.40e+6 

3 . 61e+5 

6 

7 

8 

9 

10 

11 

12 

4.76e+6 H 

3 . 60e+6 

3 . 34e+2 



4.02e+l 

2 . 05e+l 


Table 9: Condition number profile of the matrix in Example 4.3 
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classical 

look-ahead 

Example 4.1 
Example 4.2 
Example 4.3 

0.50112484863612 
0.66589668180614 
3 . 20421263631 1192e-10 

5 . 232908767834996e-15 
4 . 0288605 12358659e-14 
7 . 089249323695304e-14 


Table 10: Comparison of the relative errors for the three matrices in Example 4 


m 

0 

1 

2 

3 

4 

5 

6 

«(T m ) 

1 . 09e+15 

1.00 

3 . 98e+14 

6 . 53e+14 

8.74e+14 

1.81e+l 

7.21 


Table 11: Condition number profile of the matrix in Example 5 


look-ahead algorithm for solving general Toeplitz systems. This procedure is an extension of 
the classical Levinson algorithm for strongly regular Toeplitz matrices. We stress that our 
look-ahead algorithm is different from other generalizations of the Levinson algorithm that 
have been proposed. Except for the procedure devised by Chan and Hansen [8], all these 
algorithms can only skip over exactly singular submatrices. The algorithm in [8] allows to 
skip over ill-conditioned submatrices; however, as we showed, the procedure can break down 
if two or more consecutive blocks of ill-conditioned submatrices occur. In contrast to these 
other proposed extensions of the Levinson algorithm, our look-ahead procedure skips over 
exactly singular, as well as ill-conditioned leading principal submatrices, and it can handle 
arbitrary consecutive blocks of ill-conditioned submatrices. We reported results of numerical 
experiments, which demonstrate that our look-ahead Levinson algorithm generates solutions 
of Toeplitz systems with ill-conditioned submatrices to nearly full accuracy. 

We remark that similar techniques can be used to derive a look-ahead algorithm for 
solving general Hankel systems. A detailed description of such a look-ahead Hankel solver 
can be found in [15]. We stress that the Hankel case is fundamentally different from the 
Toeplitz case. The Hankel case is actually simpler in the sense that bilinear forms associated 
with Hankel matrices are always symmetric, and consequently one only has to deal with one 
sequence of formally orthogonal polynomials (FOPs), rather than two sequences of FBOPs 
as in the Toeplitz case. We note that FOPs associated with Hankel matrices are intimately 
connected with the nonsymmetric Lanczos process [34] for matrix computations; see, e.g., [14] 
and the references given there. 

In future work, we intend to give a rigorous stability analysis of the look-ahead Levinson 
algorithm proposed in this paper. Furthermore, we plan to develop a software package with 
FORTRAN and MATLAB implementations of this algorithm, as well as the Hankel solver 
described in [15]. 

In recent years, various so-called super-fast algorithms were developed that solve Toeplitz 
systems with only 0(nlog 2 n) operations, see, e.g., [1, 4, 5, 45]. Bunch [7] discussed the 
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stability of some of these algorithms, and he pointed out that they are unstable in the case 
of general Toeplitz systems. It is natural to ask whether look-ahead techniques can also be 
used to enhance the stability of super-fast Toeplitz solvers. This question is addressed in a 
recent paper by Gutknecht [21] who devised two super-fast Toeplitz algorithms with look- 
ahead. However, no numerical results are given in [21], and it remains to be seen whether 
these algorithms can be turned into practical numerical procedures. 
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